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Abstract 

In  navigation  satellite  systems,  it  is  necessary  to  determine  the  difference  between  the  on-board 
time  and  the  reference  time  for  each  satellite.  This  offset  can  be  estimated  in  real  time  by  filtering 
time  measurements  collected  over  a  ground  station  network,  and  has  to  be  extrapolated  when  the 
satellite  is  out  of  visibility  of  this  network.  This  analysis  has  been  carried  out  at  the  CNES  in  the 
GNSS2  and  GALILEO  context  and  leads  to  specifications  on  adjustment  and  extrapolation  errors 
of  the  on-board  time. 

The  purpose  of  this  paper  is  the  estimation  of  the  difference  between  the  extrapolation  and  the 
real  on-board  time. 

After  the  description  of  the  method,  an  example  is  given  using  real  data,  and  the  predicted 
extrapolation  uncertainties  are  compared  to  the  real  extrapolation  errors. 


INTRODUCTION 

In  1998,  CNES  proposed  a  new  concept  of  orbit  determination  and  synchronization 
for  navigation  satellite  systems  M.  In  this  architecture,  an  on-board  filter  processes 
measurements  collected  over  a  dedicated  ground  station  network  and  uploaded  in  quasi 
real-time,  in  order  to  compute  the  ephemeris  and  the  on-board  synchronization.  When 
the  satellite  is  out  of  visibility  from  the  ground  station  network,  this  synchronization 
has  to  be  predicted.  The  maximum  duration  of  extrapolation  is  3.5  hours  with  the 
station  network  considered  by  CNES. 

In  order  to  estimate  the  time  error  of  the  on-board  oscillator  when  the  satellite  is 
hidden,  a  parabolic  fit  is  performed  over  the  sequence  of  observed  time  error  data,  and 
extrapolated  during  the  hidden  sequence. 
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The  needs  of  orbit  determination  and  synchronization  were  specified  as  the  maximum 
deviation  of  the  time  error  from  the  extrapolated  parabola.  The  question  is  then:  how 
is  this  maximum  deviation  related  to  the  noise  levels  of  the  oscillator?  Several  papers 
have  already  addressed  this  question  121  131  M,  but  a  new  approach  was  chosen  here 
because  we  are  not  only  interested  in  the  asymptotic  trend  of  this  maximum  deviation, 
but  also  in  its  evolution  close  to  the  interpolated  sequence. 

In  this  paper,  we  will  denote  by  “Time  Interval  Error”  (TIE)  the  difference  between 
the  extrapolated  parabola  and  the  real  time  error  x(t)  (seeFigure  1,  left).  By  definition, 
the  TIE  samples  are  then  the  residuals  to  this  parabola' (seeFigure  1,  right). 

The  TIE  is  due  to  two  effects:  the  error  of  determination  of  the  parabola  parameters 
and  the  error  due  the  noise  of  the  oscillator.  Obviously,  both  of  these  errors  may  be 
positive  or  negative,  and  the  ensemble  average  of  the  TIE  is  equal  to  zero  (seeFigure 
2,  left).  Moreover,  it  can  be  easily  shown  that  the  ensemble  statistics  of  the  TIEare 
Gaussian  (see  figure  2,  right).  Consequently,  we  only  have  to  estimate  the  variance  of 
the  TIE  in  order  to  completely  define  its  statistical  characteristics. 

Moreover,  the  removal  of  a  quadratic  fit  from  the  time  error  sequence  cancels  out  the 
non-stationarity  problem  of  very  low  frequency  noises  (see  the  moment  condition  in 
I5!),  and  the  variance  of  the  TIE  (i.e.  the  “true  variance”)  converges  for  all  types  of 
noise  without  considering  a  hypothetic  low  cut-off  frequency. 

In  order  to  determine  an  estimation  of  the  TIE,  we  will  first  redefine  the  interpolation 
method.  Then  we  will  compare  the  equations  giving  the  theoretical  estimates  of  the 
variance  of  the  TIE  to  simulations  and  to  real  data. 


INTERPOLATION  METHOD 

Interpolating  Functions 

Let  us  consider  a  sequence  of  N  time  error  data  x(t),  regularly  spaced  with  a  sampling 
period  r0:  {x(t0),  x(<i), . . . ,  x(t;v-i)},  and  U  =  ir0. 

Rather  than  carrying  out  a  classical  quadratic  least  squares  interpolation: 

x(t)  =  Co  +  C\t  +  Cot~  +  e(t)  (1) 

where  e(t)  is  the  noise,  i.e.  the  purely  random  behavior  of  x(t),  we  use  the  first  three 
Tchebytchev  polynomials  [51  M  as  interpolating  functions  (seeFigure  3): 
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The  interpolation  we  use  is  then: 


(2) 


x(t)  —  Po$o{t)  +  +  e(f). 


(3) 
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where  the  parameters  {P0,  Pi,  P2)  have  the  same  dimension  as  x{t),  i.e.  time. 

The  classical  parameters  {Co,Ci,C2}  of  (1)  may  be  easily  deduced  from  the  parameters 
{P0.P1.-P2}  of  (3). 

Besides  their  dimensionless  nature,  the  advantage  of  using  these  interpolating  functions 
stems  from  their  normality  and  orthogonality,  which  greatly  simplify  the  estimation  of 
the  parameters  Po,  Pi, and  P2> as  well  as  their  statistical  characteristics,  as  will  be  shown. 
Moreover,  the  Tchebytchev  polynomials  minimize  the  truncation  errors,  because  their 
covariance  matrix  is  optimized  for  avoiding  bad  conditioning  problems. 


Properties  of  the  Interpolating  Functions 

Let  us  define  the  vector  associated  to  the  interpolating  function  as: 

*,-(*o)  \ 


*.■  = 


(4) 


It  is  then  possible  to  build  the  matrix  [<£]: 

/  ^o(^o)  $i(*o) 


[<J>]  =  (So  <&1  $2)  = 


$2(M  ^ 

\  $o(*jv-i)  $i(*n-i)  1)  J 


(5) 


One  of  the  main  properties  of  the  Tchebytchev  polynomials  lies  in  the  orthonormality 
of  the  vectors  associated  to  these  interpolating  functions: 


[*]T[*l=[/3l 

where  [73]  is  the  unit  matrix  (3  x  3). 

Estimation  of  the  Parameters  {P0,P1,P2} 

x{tQ)  \ 


(6) 


Let  us  define  the  vector  X  as: 


X  = 


x(tN- 1)  / 


This  vector  may  be  modeled  by: 


X  =  [$]P  +  E 


(8) 


where  P  is  the  vector  whose  components  are  the  three  parameters  we  want  to  estimate, 
and  E  the  vector  containing  the  purely  random  part  of  X: 


(9) 


Po  \ 

(  e{t0) 

P1 

and 

E=  : 

P2 ) 

\  e{tN- 1) 

In  order  to  estimate  P,  we  have  to  calculate: 

[<&]TX  =  [$]r[<I']P  +  [$]T£:. 


(10) 
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From  Equation  (6)  and  because  the  ensemble  average  =  0,  an  estimator  of  P  is 

P  defined  by: 

P  =  [$]TX.  (11) 

Thus,  the  estimate  Pj  of  the  parameter  Pj  is  easily  obtained  by  calculating: 

N- 1 

Pj  =*J-X=  (12) 

i=0 


ESTIMATION  OF  THE  TIME  INTERVAL  ERROR  (TIE) 

Estimation  of  the  Residuals 

From  (11),  the  residuals  may  be  defined  as  a  vector  R: 


R  =  X-  [<j>]P. 

The  variance  of  the  residuals  cr;  may  be  estimated  by: 
From  (13)  and  because  the  ensemble  average  RT[<p]P ^  =  0: 

(rt  ■  Rj  =  (*TX)  ~  (P'T^)  ' 


(13) 


(14) 


(15) 


The  scalar  product  (^XTx'j  is  N  times  the  variance  of  the  x(t)  data  and  the  scalar 


/  r  -\ 


product  (^PT PJ  is  the  sum  of  the  variances  of  each  estimate  P0,  Pi,  and  P2.  Thus,  the 
variance  of  the  residuals  may  be  estimated  by: 


v2e  =<r2x-J  (°Po  +  +  °Pi)  ■ 


(16) 


Correlation  of  the  Samples 

Obviously,  the  long-term  behavior  of  the  TIE  depends  greatly  on  the  type  of  noise  (from 
/~4  PM  to  white  PM).  The  autocorrelation  function  Rx{t)  of  the  x(t)  data  contains  the 
information  about  the  type  of  noise.  Since  it  is  the  Fourier  transform  of  the  spectral 
density  Sx{f).  Rx{t)  may  be  estimated  by  using: 

Rx{t)  =  2  /  cos('2nft)Sx(f)df  (17) 

Jji 

where  //  is  the  low  cut-off  frequency  and  fh  the  high  cut-off  frequency. 
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Calculation  of  the  TIE 


By  hypothesis,  we  consider  that  the  TIE  is  the  difference  between  the  true  time  error 
x (t)  at  time  t,  and  the  extrapolation  of  the  parabola  (previously  estimated  from  <0  to 
</v- 1)  up  to  this  time  t  >  t.v-o 


TIE(t)  =  x(t)  -  P0$o(t)  -  Pi$i(t)  -  P2$2(t)  and  t  >  (18) 


nru _  i-L  _ _ ] i.: ^  ~ - .hi. _ „r  a.  t'tu' - -  - i.  ~  j  i — . 

iiiua,  luc  quduidub  ciiaciiiuic  avci  age  ui  luc  jl  iu  ma/y  uc  couuuatcu  uy . 
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\*  t  y  Q  /  ^0\t)  T  y  i  /  U;  Ty  o/  ^2 1*; 

-2  [(x(t)Po)  *o(<)  +  (*(<)A)  *i(*)  +  (*(0-P2)  *2(<)] 


+2 


Y PqPi)  *o(*)*i(«)  +  / $o«)$2(t)  +  (P1P2)  $i(0$2 it) 

L  \  /  \  /  \  / 


(19) 


Consequently,  for  each  type  of  noise,  we  have  to  know:  (x2(t))  =  Rx(t),  the  autocorre¬ 
lation  function  of  x(t);  the  3  variances  (P?/  =  ap,'i  the  3  covariances  ^ PiPj ^  =  Cov(P1,PJ); 
and  the  3  covariances  (^x(t)P^  =  Cov  (x(t),  Pj). 


RESULTS  AND  DISCUSSION 

Theoretical  Results 

Since  we  are  interested  in  the  long-term  behavior  of  oscillators,  we  only  estimate  the 
TIE  for  the  3  lower  frequency  noises:  white  PM,  flicker  FM,and  random-walk  FM.  The 
theoretical  calculation  of  the  above  quantities  yields  the  following  variance: 


•  White  FM: 

•  Flicker  FM: 

(! TIE2{t )> 
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(22) 


(TIE2(t))  w  ~  0  450—  -  11102V— ^  +  93ZN2-^  -  2941V3  —  +  231V4  )  . 

315A  \  r0  r0  rj  ro  / 

All  the  above  equations  were  obtained  under  the  assumption  IV  »  1. 

Estimation  of  the  TIE  Using  the  Variance  of  the  Residuals 

The  relationships  (20)  to  (22)  need  an  explicit  knowledge  of  the  noise  levels  ka.  How- 

QTror  tfortr  1/^n  A’_f  orm  inf  ornnl  o  firm  (  o  atrot*  o.  1  /I  tvi  oti  Ka  on  /\f  f  Ka 
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type  of  noise:  the  flicker  FM  for  a  cesium  clock  or  the  random-walk  FM  for  a  quartz 
oscillator. 

Thus,  the  variance  of  the  TIE  may  be  estimated  directly  from  the  variance  of  the 
residuals  : 


•  Flicker  FM:  from  (16)  we  obtain 


•> 


7T  -N2k-3T$ 

24 


(23) 


and  from  (21)  and  (23): 
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Random-walk  FM:  from  (16)  we  obtain 


and  from  (22)  and  (25): 

2<t2 
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It  may  be  noticed  that,  for  a  fixed  value  of  <r2,  the  asymptotic  ratio  of  Equation  (26) 
over  Equation  (24)  is  exactly  equal  to  3  when  t  tends  toward  infinity. 

However,  this  method  is  less  precise  than  the  use  of  a  correct  estimation  of  the  noise 
levels,  due  to  the  statistics  of  the  estimate  of  the  variance  of  the  residuals.  We  have 
observed  experimentally  that  this  estimate  is  x2  distributed  with  a  small  number  of 
degrees  of  freedom  (3  for  a  flicker  FM  and  2  for  a  random-walk  FM).  Consequently, 
the  standard  deviation  of  the  TIE,  i.e.  the  square  root  of  equations  (24)  and  (26), 
follows  a  Student  law.  Thus,  the  bounds  given  by  the  square  root  of  Equations  (24) 
and  (26)  still  contain  68%  of  the  realizations,  but  the  95%  confidence  interval  (2a)  is 
more  than  3  times  higher. 


Comparison  with  Simulations 

In  order  to  verify  the  Equations  (20)  to  (22),  we  simulated  time  error  sequences  of 
different  types  of  noise  (white  FM,  flicker  FM  and  random-walk  FM).  For  each  type 
of  noise,  10,000  realizations  were  calculated  with  the  same  noise  level1  (fc_ 2  =  1.4  10_4s, 

3  =  3.3  10-8  or  /:_4  =  5.0 - 10~  12s— 1 ) ,  the  same  number  of  data  (65,536),  the  same  number 
of  data  taken  into  account  for  the  fit  (8640),  and  the  same  time  of  estimation  of  the 
TIE  (£  equals  to  {8639,  9900,  11350,  13000,  14900,  17000,  19500,  22400,  25700,  29400,  33700,  38600, 
44300,  50700,  58100,  65535}).  The  noise  levels  were  chosen  such  that  the  variance  of  the 
residuals  equals  1. 

‘We  denote  ka  the  noise  levels  of  the  time  error  spectral  density  Sx(f),  and  ha+2  the  corresponding  noise  levels 
of  the  frequency  deviation  spectral  density  Sy(f).  These  coefficients  are  related  by:  ha+2  =  4K2ka. 
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Figure  4  shows  the  curves  corresponding  to  the  square  root  of  equations  (20)  to  (22) 
compared  to  the  standard  deviation  estimated  from  the  10,000  simulations  (i.e.  half 
the  width  of  the  Gaussian  ofFigure  2).  The  simulations  exhibit  a  quite  good  agreement 
with  the  theoretical  curves. 

Application  to  Real  Oscillators 

Figures  5  to  8  compare  the  long-term  behavior  of  2  real  quartz  oscillators  to  the  bounds 
given  by  the  estimated  standard  deviation  of  the  TIE  (the  square  root  of  equations 
(20)  to  (22)). 

The  data  from  the  oscillators  are  time  errors  sampled  with  a  sampling  period  r0  =  10  s. 
The  fit  was  carried  out  over  the  first  24  hours  of  each  sequence  and  extrapolated  over 
the  whole  sequence  (6  days  for  oscillator  1  and  90  hours  for  oscillator  2). 

For  oscillator  1.  the  Allan  variance  yields  /i_2  =  4  ■  10-29  s-1,  i.e.  fc_4  =  lO-30  s_1.  Thus, 
the  bounds  ofFigure  5  were  obtained  by  using  the  square  root  of  (22).  On  the  other 
hand,  the  variance  of  the  residuals  of  the  interpolated  part  (i.e.  the  first  24  hours) 
is  =  1,4-  10-17  s2.  Thus,  the  bounds  ofFigure  6  were  obtained  by  assuming  that  the 
random-walk  FM  was  dominant  and  by  using  the  square  root  of  (26). 

For  oscillator  2,  the  Allan  variance  revealed  that  2  types  of  noise  must  be  taken  into 
account:  white  FM  (h0  =  10~22  s,  i.e.  A-_2  =  2,5-  10-24  s)  and  random  walk  FM  (/i_2  = 
1, 5  -  10— 31  s_1,  i.e.  &_4  =  4  •  10-33  s-1).  Thus,  the  bounds  ofFigure  7  were  obtained  by 
using  the  square  root  of  the  sum  of  (20)  and  (22).  On  the  other  hand,  the  variance  of 
the  residuals  of  the  interpolated  part  (i.e.  the  first  24  hours)  is  <r2  =  1,4-  10— 18  s2.  In 
this  case  also,  we  assumed  that  the  random-walk  FM  was  dominant,  and  the  bounds 
ofFigure  8  were  obtained  by  using  the  square  root  of  (26). 

The  experimental  TIE  curves  remain  in  the  theoretical  bounds  except  for  Figure  6. 
This  is  quite  compatible  with  the  statistics  of  TIE  since  32%  of  the  realizations  should 
be  outside  the  bounds. 


CONCLUSION:  A  New  Strategy  for  Long-Term  Stability  Analysis 

Besides  the  interest  of  this  method  for  synchronization  prediction,  it  may  also  be  used 
for  defining  a  new  method  for  very  long-term  stability  analysis. 

An  oscillator  may  be  continuously  measured  during  a  few  days  (e.g.  a  time  error 
measurement  with  a  sampling  period  of  1  minute  during  10  days).  From  these  data, 
the  noise  levels  of  this  oscillator  could  be  precisely  determined  t7)  t6l  and  a  quadratic 
fit  could  be  carried  out.  Thus,  if  the  oscillator  is  continuously  running  in  the  same 
conditions,  it  could  be  possible  to  extrapolate  the  difference  of  this  oscillator  with  the 
parabolic  fit  after  a  few  months  or  one  year. 

This  analysis  could  be  helpful  in  low  accuracy  time  keeping  applications,  for  instance 
for  industrialists  who  periodically  send  their  oscillator  to  an  accreditation  laboratory. 
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Figure  1:  Quadratic  fit  over  a  time  error  sequence  of  white  FM  (left).  The  fit  is  performed  over  the  first 
quarter  of  the  sequence.  The  residuals  of  this  fit  (right)  correspond  to  our  definition  of  the  TIE. 


Figure  2:  Dispersion  of  the  TIE  for  20  realizations  of  the  same  white  FM  process  (left).  The  histogram  of 
10000  realizations  of  the  TIE  estimates  (here  for  t  —  40,000s)  exhibits  a  Gaussian  behavior. 


'IA'1 


J1J 


0.25 


Figure  3:  The  first  3  Tchebytchev  polynomials  calculated  for  N  =  100  data 


Figure  4:  Comparison  of  the  estimation  of  the  standard  deviation  of  the  TIE  calculated  from  the  equation 
(20)  to  (22)  (solid  lines)  and  estimated  over  10,000  realizations  of  simulated  noise  (circles,  squares,  and 
triangles).  The  lower  curve  was  obtained  with  white  FM  (circles),  the  middle  one  with  flicker  FM  (squares), 
and  the  upper  one  with  random  walk  FM  (triangles).  In  order  to  use  the  same  scale,  the  noise  levels  were 
defined  in  such  a  way  that  the  variance  of  the  residuals  is  equal  to  one  (fc-2  =  1-4  •  10~4s,  k-3  =  3.3  ■  10-8, 
i t_4  =  5.0  •  10-12s-1).  The  error  bars  corresponding  to  the  estimates  of  the  simulated  noises  are  too  small 
to  be  plotted  on  this  graph. 
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Figure  5:  Oscillator  1.  Evolution  of  the  TIE  after  the  fitted  sequence.  The  fit  was  performed  over  1  day  with 
a  sampling  period  of  10s.  The  estimation  of  the  TIE  (dashed  line)  was  performed  from  Equation  (22)  by 
using  a  noise  level  estimate:  Ar_4  =  10-30s-1.  On  the  right,  the  plot  is  an  enlargement  of  the  first  4  hours. 
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Figure  7:  Oscillator  2.  Evolution  of  the  TIE  after  the  fitted  sequence.  The  fit  was  performed  over  1  day 
with  a  sampling  period  of  10s.  The  estimation  of  the  TIE  (dashed  line)  was  performed  from  Equations  (20) 
and  (22)  by  using  the  noise  level  estimates:  fc_2  =  2,5  •  10-24  s  and  fr_4  =  4  •  10~33  s-i.  On  the  right,  the 
plot  is  an  enlargement  of  the  first  4  hours. 


Figure  8:  Oscillator  2.  Same  as  fig.  7.  The  estimation  of  the  TIE  (dashed  line)  was  performed  from 
Equation  (26)  by  using  the  variance  estimate  of  the  residuals  cr2  =  1,4  •  10~18  s2.  On  the  right,  the  plot  is 
an  enlargement  of  the  first  4  hours. 


